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Abstract 

Thermal oscillations of base pairs in the Peyrard-Bishop-Holstein model are 
simulated by stochastic fluctuations of base overlap integrals. Numerical 
investigation of the model is carried out for a hole transfer in G1A2G5G4 se- 
quence which was previously studied experimentally by F. Lewis et al. A hole 
migration between Gi and G3G4 is determined by the matrix elements of the 
charge transition, but presence and amplitude of their stochastic fluctuations 
proved to play a key role in reproduction of the experimental kinetics. Good 
agreement with the experimental data was obtained for a wide range of the 
model parameters' combinations. 
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1. Introduction 

Investigations of DNA conducting properties are very important for both 
classical radiobiology [l], ^ and quite a new science of nanobioelectronics 
js], 0] . Charge migration in DNA is strongly dependent on a set of conditions, 
the basic of which is the nucleotide sequence. Researchers are generally 
focused on a cation-radical (hole) migration. A detailed consideration of 
experimental and theoretical material on charge transfer in DNA is presented, 
for example, in reviews by Conwell jsj and Wagenknecht j6[. 
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Experimental investigations of a hole transfer in DNA can be divided into 
steady-state and time-resolved measurements. The first ones are based on 
ionization of modified DNA followed by the analysis of its oxidation prod- 
ucts by electrophoresis, HPLC or other methods. Owing to steady-state 
methods it has been found, that a cation-radical can move in DNA over a 
distance of tens or even hundreds of angstroms 0] ^ @] • Theoretical studies 
demonstrated that in the case of heterogeneous DNA chains, a hole transfer 
is realized as a series of hops between guanine nucleobases, which have the 



lowest oxidation potential [10| - [12|. Theoreticians discuss two basic models 



of the cation-radical transition, which were originally suggested by J. Jortner 



et al. [13|. In the first, superexchange model, a charge tunnels immediately 
from a donor to an acceptor, avoiding chemical interaction with interven- 
ing nucleotides. The second, multistep hopping model, implies a successive 
transition of a hole from one base to another due to the energy of thermal 
fluctuations. In the case of regular and homogeneous chains the transition 
can also be realized by a band, polaron, soliton, breather, etc. mechanisms 
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Though the equilibrium methods have played a major role in the study 
of charge transfer, the information of relative hopping rates is incomplete 
for estimating the time of a hole migration over a particular distance. By 
contrast, time-resolved measurements enable one to get absolute rates of 



separate charge hopping steps in a DNA fragment [15|] - [18 



Prominent among all time-resolved measurements are the works by F. 



Lewis et al. ly, ll7|. Unlike many other works (see, for example, 15|, ll8|), 
there the object of investigation was oligonucleotides shorter than 10 base 
pairs. One end of these oligomers was bound to stilbene-4,4'- dicarboxam- 
ide chromophore (hereafter St) which actually made them hairpins. A hole 
transfer in such structures takes place in nanosecond timescale. The general 
form of hairpins from work [l6j was St-Am-GiA2G3G4,-An, where 1 < m < 3, 
and < n < 2. Hereafter, 1-4 are numbers of bases in the "key fragment". 
The singlet stilbene selectively photooxidizes Gi, resulting in the formation 
of a primary radical-ion pair {St~* Gf*A2G3G4}. In the subsequent 
reversible transfer of a hole to the G3G4 doublet Gf* in turn acts as an elec- 
tron acceptor. Hence, both a donor and an acceptor are guanine bases in 
the hairpins under study. This fact together with a small size and simple 
primary structure of the hairpins leads to a simple reaction scheme. The 
scheme enables direct assessment of the rate constants for a hole transfer 
from Gi to G^G^ and back. 
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Due to a simple and elegant technique F. Lewis et al. were the first 
to determine in 2000 the values for rate constants of the forward and return 
cation-radical transport {kt and k-t respectively). Their values were 50 and 6 
lis~^ respectively. These data are in g ood agreement with other time-resolved 



measurements (see, for example, 18| ) 



Nevertheless, the theoretical interpretation of the charge transfer dynam- 
ics presented in jToj is based on a phenomenological description with the use 
of a reaction scheme. It does not make possible studying physical principles 
of the cation-radical migration in heteropolymer DNA. The aim of this work 
is to describe the results of 16| by a microscopic model and explain them 
relying on the data obtained in the computational experiment. 

Here we studied numerically the Peyrard-Bishop-Holstein model for a 
fragment of a G1A2G3G4, sequence. The model allows an accurate taking 
account of the base-pair dynamics in terms of the Peyrard-Bishop-Dauxois 
Hamiltonian, which has been sucessfully applied for the study of various 



aspects of DNA denaturation for over 20 years [19| - [21|. In our approach 



thermal oscillations of the lattice structure are introduced via small stochastic 
fluctuations of nondiagonal matrix elements of the charge transition in the 
quantum subsystem. 

These matrix elements present the only group of parameters that was 
obtained in quantum-chemical computations (see refs below). This fact is 
very important: taking into account the influence of the temperature, hy- 
dration and short length of the duplex, nondiagonal elements can possess 



some other values in real DNA in water solution 22, 231. Hereafter the term 



"parameter combination" implies combination of values of the nondiagonal 
matrix elements and amplitudes of their small stochastic fluctuations. All the 
other parameters of the Peyrard-Bishop-Holstein model are experimentally 
measured quantities and remained invariant in our simulations. 

In the numerical experiment, a hole transfer over the G1J42G3G4 fragment 
is investigated for more than 100 parameter combinations. These simulations 
enabled us to find a set of parameter combinations at which the charge trans- 
fer rates averaged over the " microensemble" are in good agreement with the 
experimental data obtained by Lewis and co-workers. So, the experimental 
data on the cation-radical transfer along DNA are for the first time repro- 
duced by microscopic modeling. Below we describe in detail the model and 
the approach developed, discuss our results and compare them with the ex- 
periment. 
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2. The model 



The Peyr ar d- Bishop- Holstein (PBH) model is a recently developed "hy- 



brid" of the Holstein [2J] and Peyrard-Bishop-Dauxois [19[ models. So far, 
the PBH model has been used to study the charge migration only in ho- 
mopolymer DNA 25|] - 27|. In the studied case of a heterogeneous sequence, 



the motion equations for a hole and nucleotide pairs in the neighborhood ap- 
proximation have the form: 



n-1 + ^n,n+l 



dt 



OUn dt 



(1) 



In the first equation, h is the Planck constant, ipn — the probability am- 
plitude for the charge carrier located at the n-th site, — its oxidation 
potential, S' is the charge- vibrational coupling constant, represents the 
transverse displacement of the hydrogen bonds connecting two bases, i'n,n±i 
are matrix elements of the transition. In the second equation, m is the effec- 
tive site mass, 7 — friction constant and 



ViUr, 



DJe 



(2) 



W{un-i^un) = ^(l + pe-'^(""+^"-))K-S„_i)2 (3) 

where D„ and determine the depth and width of the energy well in Morse 
potential V{un)- In the nearest-neighbor potential W{un-i,Un), accounting 
for stacking interactions, k — coupling constant, /3 — damping coefficient 
and p is dimensionless stiffness parameter. The condition p 7^ reproduces 
the fact, that the double-stranded backbone is more rigid than the unwound 
strands. 

The friction constant was taken to be equal to 6 ■ 10"^"^^^ • rrT^ ■ s, which 
corresponds to the picosecond characteristic time scale of DNA oscillations 



28l |. The other parameters of the classical subsystem were taken from the 



well-known Letter by Campa and Giansanti: -D^r = 0.05 eV, Dec = 0.075 
eV, Oat = 4.2 A-\ ace = 6.9 A'^, k = 0.025 eV ■ A'^, /3 = 0.35 A-\ p = 
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2 29| . The oxidation potentials for adenine and guanine were taken from 
experimental data for acetonitrile solution: v\ = 1.69 eV, Uq = 1.24 eV 30 



The potential for G, being the lowest one, was taken as zero and, hence, i^^ 
= 1.69 — 1.24 = 0.45 eV. The estimated values of the matrix elements on 
the G1A2GSG4 fragment were also taken from literature: uga = 0.089 eV, 



z/^^ = 0.049 eV and i^gg = 0.084 eV [22 



When transforming system (1) into dimensionless form, the condition was 
specified: 



m ■ U h 

where r — the timescale of the system, U is an arbitrary scale of the trans- 
verse displacement , x ^ the dimensionless form of the charge- vibrational cou- 



pling constant. If we choose 6' = 0.13 eV ■ 28|, then U = 10~^^ m and 
X ~ 0.02, provided that m = 10~^^ kg and r = 10~^^ s. Consequently, the 
parameters of the model relate to their dimensionless forms as 

d„ = an ■ U~^] t = t-T] P = /3 ■ U~^] Un = Un-U] 

Dn = Dn ■ m ■ U"^ ■ r^^; j = u' ■ m ■ t^-^; k = Vl? ■ m ■ r~^; 

where ?7n,n±i) Vni ^' ^ — dimensionless values of Vn^n±\-, 7 ^^"^ ^ 
respectively. 

The dimensionless form of system (1) was realized in a parallel MPI- 
program. Numerical calculations were carried out by the fourth-order Runge- 
Kutta method, the step size being 10~^^ s. Each separate realization was 
simulated for over 30 ns (3 ■ 10^° steps). At the initial moment the charge 
was on Gi, -01 = 1. To specify boundary conditions we introduced fictitious 
terminal base pairs Xq and X5, for which i'n,n±i and Un had a fixed value — 
zero. Hence, hole migration to Xq and X5 in the X0G1A2G3G4X5 sequence 
was impossible. Such boundary conditions were referred to as fixed. 

3. The method 

In a real duplex, thermal vibrations of the bases lead to fluctuations of 
overlap integrals of their vr-orbitals. Since DNA is a complex molecule with 
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a generous amount of vibrational modes, the time dependence of ?7n,n±i is a 
very complicated stochastic function. Nevertheless, at this stage of model- 
ing a key moment is taking account of fluctuations of non-diagonal matrix 
elements itself. Actually, a more careful consideration of these perturbations 
is desirable, but in the first approximation it is enough to specify the only 
(constant) frequency of ?7„,n±i fluctuations. We take this parameter to be 
equal to the frequency of the stretched oscillations (along H-bonds), i. e. 
0.5 • 10^^ s~^. In this case, the amplitude of each individual oscillation is a 
random value, whereas the root-mean-square amplitude of the fluctuations 
is a parameter of the model. 

In our numerical experiments thermal oscillations were simulated by ran- 
dom deviations of the elements rjcA, Vag and tjgg from their "reference" 
values specified as parameters. The reference value of any matrix element 
(let us call it r]i, ) was invariable throughout the realization. Fluctuations 
were provided by the stochastic addition 77. In such a way, the resulting 
value of the matrix element was the sum r/^ + r]. Generated pairs of random 
numbers ri and r2 whose values were in the range from to 1, were then 



subjected to Box-MuUer transformation [SjJ- This resulted in 77^ — a random 
addition r] in the next point of extremum 



Va 



S ■ cos(27rri) ■ a/— 2 In r2 



where S is the root-mean-square amplitude of a random deviation. 

At initial moment r] was equal to zero. At that very moment rja was 
generated for the next fluctuation: rja = ri{t + 10~^^s). In the course of the 
simulation r] changed with the step j = 10"^^ s according to the law 



V{t + J) = — ^ + 2 cos{n- J -10 s ), 

where 7]'^ is the value of the previous deviation (initially equal to zero), and 
j = 1 ■ 10-1^ s, 2 ■ 10-1^ s, . . . , 99 ■ 10^1^ s, 10^12 At the moment j = IQ-^^ 
s r]'^ reached a value of rja, while rja was generated anew. 

For any combination of the matrix elements, the values of their standard 
deviations So a, Sag and Sgg always related to one another by the constant 
ratio of 7 : 5 : 8. This was caused by the ratio of the quantities rjcA, Vag 
and rjGG per se (8.5 : 4.7 : 8). A slightly reduced value of Sga reflects the 
fact that Gi resides in the center of the duplex and is subjected to strong 
fluctuations less than the other base pairs under consideration. For each 
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combination of the parameters Sqa^ Sag^ Sggj Vga, Vag and tjcg we per- 
formed 80 realizations to obtain the curves ( I'j/'iP) = fU), where angular 

\ /so 

brackets denote averaging. 
4. Results and discussion 



According to the data by A. Voityuk et al. |22[, tjga = 1-351, tjag = 
0.743, and tjgg = 1-275. Taking into account hydration of oligonucleotide, 
we specified variation of these quantities to be approximately 25%. Thus, 
the reference values (?7b, see above) of the matrix elements was: r]GA = 
1.051,1.351,1.651; r]AG = 0.543,0.743,0.943; r]GG = 1-275,0-975- We ex- 
cluded variant rjcG = 1-575 because of so-called "end-fraying", causing sub- 



stantial reduction of corresponding overlap integral |32| - [3J] - The maximal 
standard deviations of the nondiagonal matrix elements was chosen such, 
that the probability of a situation rj^ + rja < he less than 10~^- 

In the first series of simulations we took maximum dispersion of the matrix 
elements of the charge transition: Sqa = 0-14, Sag = 0-10, Sqg = 0-16. The 
transfer was studied during 14 ns. For r]GG = 1-275, in all the cases, most 
of the charge density remained on Gi throughout the realization. The only 
exception was combinations of the parameters tjga = 1-651, tjag = 0-743 and 
VGA = 1-651, rjAG = 0.943, where very slow transfer was observed. In the case 
of Vgg = 1-025 the charge escape from the first site took place for all the 
chosen rjcA and rjAG though its rate was quite low. We hypothesized that for 
any combination of rjcA and r]AG (in some range of these parameters) there 
exist a tjgg and Sga, Sag, Sgg such that the rates of direct and backward 
transfer of a hole over the G1742G3G4 fragment are close to the experimental 



values obtained by Lewis et al, i.e. kt ~ 50 fis~ , k_t ~ 6//s~ [16 . 

In the second series of simulations we extended the time up to 30 ns and 
studied the model using 41 combinations of Vga, Vag and ijgg for strong 
fluctuations {Sga = 0.14, Sag = 0.10, Sgg = 0.16), using 54 combinations 
for moderate fluctuations {Sga = 0.07, Sag = 0.05, Sgg = 0.08) and using 27 
combinations for weak fluctuations {Sga = 0.035, Sag = 0.025, Sgg = 0.04). 
The solution of kinetic equations, describing the process 

G'i"*^2G^3G^4 ^ GiA2G3G^* (4) 

which had the form 
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was used to fit the curves ( "01 / = fit) obtained for each combination 

\ /so 

of the parameters. Constants kt, k^t and the constant of proportionahty F 
were chosen such that the root-mean-square difference between approx 

and ( IV^iP) (t) on the time interval from 10 to 30 ns be minimum. 

\ /so 

For one of the functions ( ) (t), together with finding the constants 

\ /so 

kt and k_t, we carried out comparison with experimental data by Lewis. 
In the kinetic scheme of Letter [16| consideration is given to two states of 
an ionized hairpin, namely, a radical-ion pair X having the form of {St~' 
G^' 7426*3 6*4} and capable of recombination and a radical-ion pair Y {St~' 
GiA2G3G^*}, which cannot recombine. In our computational experiment a 
recombination reaction is lacking, therefore a straightforward comparison of 
our results with the experimental data by Lewis et al. is impossible. Never- 
theless, knowing [X]{t) / {[X]{t) + Y{t)) = (t) at every instant, we 
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can carry out "indirect" substitution of our data into the kinetic scheme from 



Letter [16[. Comparison of the computational data with the experiment is 
presented in Fig. 1. 

For almost all combinations of "qcA, Vag and tigg: the closest to Lewis' 
values of kt, k_t and equilibrium constants were obtained for intermediate 
standard deviations: Sga = 0.07, Sag = 0.05 and Sgg = 0.08. These 
data are presented in Table 1. It turns out that S plays a key role in the 
relation between kt and k-t, obtained in the calculations: strong fiuctuations 
of the matrix elements lead to a decrease in kt and increase in k_t, that is 
to a decrease in the equilibrium constant. For example, the combination 
of the parameters Sga = 0.14, Sag = 0.10, Sgg = 0.16, tjga = 1-351, 
Vag = 0.743, tjgg = 0.705, yields kt = 29fxs~^, k^t = 20fis"^: the equilibrium 
constant turns out to be almost 6 times less than the experimental value. 
Similarly, too weak fiuctuations lead to abnormally large (more than 20) 
equilibrium constants due to negligible k_t: 0.01 - 0.5 yus~^. Moreover, a 

standard deviation determines the dynamics of ^IV'iP^ (^) decay during 
the first few nanoseconds. At large S the initial escape of a hole from Gi 
occurs very fast [kt > 100/is~^ during the first 3-4 ns). In the case of too 
weak fiuctuations, by contrast, the initial rate of a hole migration is rather 
low as compared to its subsequent transfer. 

The data of Table 1 suggest that one cannot unambiguously determine 
the values of the matrix elements of the charge transition at which kt and 
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Figure 1: Comparison of computational data obtained for the parameter set Sga — 0.07, 
Sag = 0.05, Sgg = 0.08, tjga = 1-051, t]ag = 0.543, 7700 = 0.519, with the experimental 
data found in [l6|. The estimated values of the constants: kt = 60/is~^, = 
[lij . The best fitting to Lewis' data was achieved when kcr = 470/is~^, the quantum 
yield and recombination being taken into account too. The solid line (gray) represents the 
experimental decay curve of St~' from [lij . where 6 A is a relative transient absorption at 
575 nm. The data of the computational experiment are presented as the dot set (black), 
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Time interval between the dots is equal to 2 ps. 
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Table 1: The values of rate constants kt and fc_t, obtained for various combinations of tjga, 
TjAG and r]GG. For each combination of tjga, Vag, the value of rjGG is presented such that 
kt and (which are also shown) are closest to the experimental data [l^. " For the 
rjGA = 1-651, rjji^G = 0.943, r]GG = 0-9: a great deviation of k_t from the experimental one 
results from insufficiently high amplitude of ry fluctuations: for Sga = 0.14, Sag = 0.10 
and Sgg — 0.16 the values of kt and ar e eq ual to 32 and 8.6 ^jLs^^, respectively, which 
is slightly closer to the data by Lewis et al [16| . 



Vag 


0.543 


0.743 


0.943 


VGA 










1.051 


rjGG = 0.52 ± 0.02 


Vgg 


= 0.56 ±0.02 


r]GG = 0.56 ± 0.02 




kt = 59.2/xs-i 


kt = 


53.9/xs-i 


kt = 67.5/xs~-'^ 




k^t = 7.2//s-^ 


k-t-- 


= 6fis~^ 


k^t = 7fis~^ 


1.351 


r]GG = 0.675 ± 


Vgg 


= 0.71 ±0.02 


r]GG = 0.73 ± 0.02 




0.015 


kt = 


54.5/xs~^ 


kt = 48.3//s-^ 




kt = 52.2//S-1 


k-t-- 


= 5.5/is~^ 


k-t = 7.8/is-^ 




k^t = 6/is~^ 








1.651 


r]GG = 0.825±0.02 


Vgg 


= 0.87 ±0.02 


Tj^^ = 0.9 ±0.03 




kt = 52.5iJ,s~^ 


kt = 


46.4//S-1 


kt = A2fis-^ 




k_t = 5.85/is~^ 


k^t-- 


= 4/is~^ 


k^t = 0.15//S-1 
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k^t fits the experimental data. Moreover, the range of ?7gA; Vagj and rjcG 
combinations for which the rate of the charge migration corresponds to the 
experiment is, probably, rather wide. 

It may seem strange that all the values of ijgg in Table 1 are 1.5 - 2.5 
times less than the theoretical one (1.275 dimensionless units). However 
this deviation can easily be explained by end effects. In this respect very 
illustrative are NMR- investigations of local "flipping-out" of bases in DNA 
oligomers. Due to the end-fraying, the rates of flipping of terminal base pairs 



are so high that they cannot be measured [32| - [3J| . Besides, the influence of 



the end- fraying can extend to the second and even the third base pair [32|, |33 



In experiments by Lewis et al. end-fraying could be responsible for reduced 
overlapping of two terminal base pairs. On the other hand, deviation of 
rjGG from the literature value may be to some extent caused by imperfection 
of our calculation technique. Supposedly, if we change the flxed boundary 
conditions, for example, to the periodical ones, we will manage to get the 



values of r]GG more close to the estimated data by Voityuk et al. [22 . 

In this Letter we studied migration of a cation-radical in the G1A2G3G4 
fragment in the framework of a proposed semistochastic variant of the PBH 
model. We have succeeded to show the existence of a wide range of the 
model parameter combinations at which the rates of direct and reverse charge 
transfer are close to the experimental data. In the studies of a hole transfer 
in the framework of a more simple microscopic model, the time averaged 
rate of a hole transfer from Gi to G^G^ was approximately two orders of 



magnitude higher than in the experiment [35[. Moreover, the shape of the 



curve iV'iP = f{t) per se reminded exponential decay only slightly. In our 

computational results, only the function ^iV^iP^ = f{t) has an exponential 

form, while IV'iP = /(O for each individual realization has nothing to do 
with an exponent. Due to features of Morse potential (see eq. (2)) a positive 
charge cannot substantially deform the molecular lattice (in this case m„ < 0). 
This dramatically constraints its mobility. Therefore, the "motor" of a hole 
migration is transient and rather rare " successful" combinations of the matrix 
elements of the charge transition. 

The problem of how the transfer of a cation-radical relates to fluctuations 
of base overlap integrals is to be studied in the future. Further development 
of our approach implies the direct coupling of matrix elements' fluctuations 
with thermal dynamics of DNA in PBH model (Langevin approach, see, e. 



36|). An alternative way is the design of a more quick semistochastic 
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algorithm of charge transport simulation in a long heterogeneous DNA. For 
example, we intend to model a hole transfer in other sequences studied by 



Lewis et al. 17 



Microscopic modeling of DNA has played a great role in the studies of 



DNA thermal and mechanical denaturation 19| - 2l|]. In the same way it is 
promising for investigation of charge migration. For example, in thermody- 
namical calculations which are usually used for interpretation of experiments 
on a hole transfer in DNA, an appropriate taking account of chaotic fluctua- 
tions of overlap integrals is hardly feasible. Nevertheless, an important role 
of these fluctuations in providin g th e charge transfer along heteropolymer 
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DNA was proved experimentally 

Here a computational experiment enabled us to demonstrate an important 
role of stochastic oscillations of the bases in the cation-radical migration 
process. Moreover, consideration of these oscillations turned out to be of 
crucial importance for reproduction of experimental data on charge transfer 
in heteropolymer DNA. 
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